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Abstract 

Path integral Monte Carlo simulations are applied to study dense atomic hydrogen 
in the regime where the protons form a Wigner crystal. The interaction of the 
protons with the degenerate electron gas is modeled by Thomas-Fermi screening, 
which leads to a Yukawa potential for the proton-proton interaction. A numerical 
technique for the derivation of the corresponding action of the paths is described. 
For a fixed density of r s = 200, the melting is analyzed using the Lindemann ratio, 
the structure factor and free energy calculations. Anharmonic effects in the crystal 
vibrations are analyzed. 
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1 Introduction 



Recent advances in high pressure experiments and computer simulations tech- 
niques have led to substantial progress in our understanding of the melting 
properties of solid hydrogen at high pressure. Using diamond anvil cell experi- 
ments, E. Gregoryanz et al. [1] extended the experimental determination of the 
melting line from 15 GPa [2] to 40 GPa. Bonev et al. combined the two-phase 
melting technique with ab initio simulations [3] and predicted that there is a 
maximum in the melting temperature around 80 GPa [4]. At these pressures, 
hydrogen is still in molecular form. The question of interest is whether the 
melting temperature keeps decreasing as the pressure is increased further, or 
if another phase appears and the melting temperature again increases. 
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From theoretical arguments, we know that atomic hydrogen eventually trans- 
forms into a Wigner crystal where the protons form a body centered cubic 
(b.c.c.) lattice. The focus of this work is to characterize this regime more ac- 
curately, to determine its melting line and to analyze by how much the known 
molecular solids phase and the Wigner crystal phase are separated on the 
pressure or density scale. 

In this article, we study the Wigner crystal of protons using path integral 
Monte Carlo (PIMC). This technique allows us to characterize the quantum 
effects of the protons. The anharmonic effects in the lattice vibrations are also 
included accurately. While at much higher temperature and lower density, one 
can describe the electrons from first principles [5]. For this work, we instead 
approximate the electron-proton interaction using Thomas-Fermi theory. This 
leads to an effective Yukawa potential for the proton-proton interaction, 

V Y (r) = -e- r / D * , (1) 
r 

where the screening length D s is given by [6], 



and r s is the Wigner-Seitz radius, ^irr^ = V/N. Throughout this work, we will 
use units of nuclear Bohr radii, a = 4:7ie h 2 / m p e 2 = 2.9xlCT 14 m and nuclear 
Hartrees, Ha = e 2 /(47re ao) = 8.0xl0~ 15 J = 5.8 x 10 8 K k b , which are by a 
factor by m p /m e = 1836 shorter, or larger respectively, than the usual atomic 
units. 

Jones and Ceperley [7] used PIMC to study quantum melting in Coulomb 
systems where the electrons were assumed to form a rigid background. We 
extend their work by introducing the Yukawa potential in order to understand 
how the electronic screening affects the stability of the Wigner crystal. 



2 Path integral Monte Carlo 



The thermodynamic properties of a many-body quantum system at finite tem- 



perature can be computed by averaging over the density matrix, p 
l/k h T. Path integral formalism [8] is based on the identity, 



e~ p6 = 



— -E-H 

e M 



M 



(3) 
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where M is a positive integer. Insertion of complete sets of states between the 
M factors leads to the usual imaginary time path integral formulation, written 
here in real space, 



p(R,R';/?) = J... J 



dRi . . . dR M -i p(R, Ri; r) . . . p(R M -i, R'; r) , (4) 



where r = /3/M is the time step, and R is a collective coordinate including all 
particles, R = {ri, . . . , r^}. Each of the M steps in the path now has a high 
temperature density matrix p(Rk, Rfc+ij t) associated with it. The integrals 
are evaluated by Monte Carlo methods. For the densities under consideration, 
we can neglect exchange effects of the protons and represent them by dis- 
tinguishable particles. Given these constraints, PIMC is an exact technique 
and free of uncontrolled approximations (assuming the Yukawa potential is 
valid). This technique includes the correct phonon excitations in the presence 
of anharmonic effects, which we will discuss later. 

2. 1 Action for an Isolated Pair of Particles 

The action plays a central role in PIMC since it determines the weights of 
paths. We will describe a novel approach for its derivation, which we found 
to be more accurate for Yukawa systems than previous techniques. First we 
discuss the action for an isolated pair of particles, and then we introduce 
periodic boundary conditions commonly used in many-body simulations. 

Typically, one approximates the high-temperature many-body density matrix, 
p(R, R'; r), as a product of exact pair density matrices which can be motivated 
using the Feynman-Kac (FK) formula, 



where Po(R> R'; r ) is the free particle density matrix. Uyi^ij, •; r) is the pair 
action corresponding to all paths separated by at imaginary time t = 
and by r^- at t — r. An approximation is introduced when one makes the as- 
sumption that the different pair interactions can be averaged by independent 
Brownian random walks, denoted by brackets (...). The pair action approx- 
imation is exact for two particle problem. However, higher-order correlations 
are left out, which must be recovered in the many-body PIMC simulations 
using a sufficiently small time step r. 




= (n e ^° Ta!rtV(rii) 

\i<j 

-V. . [/y(r^,r' .;t) 



R^R' 



(5) 



(6) 
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The pair action, Uy, can be computed by three different methods. 1) For cer- 
tain potentials where the eigenstates are known in analytical form, e.g. for 
the Coulomb potential, the action can be derived from the sum of state [9]. 
However to our knowledge, for the Yukawa potential they are not known ana- 
lytically. 2) In the matrix squaring technique [10], one represents the density 
matrix on a grid and successively lowers the temperature by performing a one- 
step path integration. This method can be applied to arbitrary potentials. 3) 
For the Yukawa potential, we found it advantageous to use the FK formula 
to derive the pair action. Computationally, it is a bit more expansive than 
matrix squaring but it does not introduce grid errors that we found difficult 
to control in case of the Yukawa potential. The FK approach is also applicable 
to arbitrary potentials unless they exhibit an attractive singularity, which is 
discussed further in [11]. 

In FK approach one derives the action Uy(r, r'; (5) stochastically by generating 
an ensemble of random paths according to the free particle action Uq{y, r'; /3) 
that begin at r and terminate at r', 

e -(U Y -U ) = /g-r^X^IMri-O+WCrOA 

\ / Uo(t,t';0) 



(7) 



The free-particle paths can be generated by a bisection algorithm [12]. We 
found that M = (3/r = 32 time steps was sufficient for the Yukawa potential. 

To determine the kinetic energy in simulations, one also needs the derivative 
of the action with respect to (3, which can be evaluated from the same set of 
paths, 



M 



A (tv _ tf0) = / e -E.^ E 



i=l 



W(r ! ) + -(r ! -r^).VW(r ! ) 



(8) 



t/o(r,r';/3) 



where rf 1 represents the classical path between the two end points. 



The FK formula yields the action for only one specific pair of r to r'. For the 
diagonal action, r = r', we map out a whole grid of points beginning from 
a small value near the origin to large values (several times the thermal de 
Broglie wavelength given by 



' 2irh 2 f3 / m) . Typically, we use a logarithmic grid 
with about 500 points. For large r, the action approaches the classical limit 
given by the primitive approximation, 



U Y (r,r>;P)*^[V Y (r) + V Y (r>)} , 



(9) 



which is shown in Fig. 1. 
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Particle separation r [a B ] 

Fig. 1. The pair action Uy (r, r' = r; r) scaled by by 1/r is compared with the corre- 
sponding Yukawa potential as a function path separation, r. The Yukawa screening 
length of D s = 387 represents the electronic screening of the proton interaction in 
hydrogen at a density of r s = 200. The action converges to the primitive approxi- 
mation, Eq. 9, for large r and small r. At small r, quantum fluctuations remove the 
singularity present in the Yukawa potential and lead to a linear dependence on r, 
which is known as the cusp condition. 

Statistical uncertainties in the resulting action are intrinsic to the FK ap- 
proach. We found that 10 6 paths yield sufficiently small error bars. However, 
any noise in tabulated action values is impractical for the subsequent inter- 
polation in PIMC simulations. To eliminate this problem, we use the same 
random paths for all grid points in the table. This does not remove the un- 
certainty but prevents noise in the tables. 

Including off-diagonal density matrix elements in PIMC simulations allows one 
to use larger time steps, which makes the simulation more efficient. However, 
the off-diagonal terms are more difficult to obtain with the FK approach, 
which is one of the limitations of the approach. Here, we only consider the 
leading term in an expansion of the action [12], 

U Y (v, r'; (5) « U Y (q, q; (3) + s 2 £(g; /?) + ... (10) 



where q = |(|r| + |r'|) and s = |r — r'|. It turns out that it is more efficient 
to derive £ from finite differences in s rather than evaluating an analytical 
expression. Having completed our derivation for the action of an isolated pair 
of particles, we now consider a system with periodic boundary conditions. 
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2.2 Pair Action in Periodic Boundary Conditions 



In our PIMC simulations, we use N = 250 particles, which are initially placed 
on the sites of a b.c.c. lattice. For the density under consideration (r s = 200), 
the corresponding screening length for hydrogen (D s = 387.667) is comparable 
in magnitude to the length of the simulation cell L = 1218.59. Consequently, 
long-range effects from periodic image particles are far from negligible, and 
significant care must be taken to derive results in the thermodynamic limit. 

The total potential energy for a system of N particles interacting via the 
Yukawa potential Vy(r) is given by [13], 

i>j L Z i L^O 



where = — rj and L is a lattice vector. Instead of representing long-range 
terms, Vp(r) = ElW(f + L), on a 3D table [14], we adopt the optimized 
Ewald [15], technique by Natoli and Ceperley [16] and express the potential 
as a sum of one real-space image, W(|r|), and a number of Fourier components, 



Vp(r)=J2V Y (r + L)^W(\r\)+ ]T y k e +l 



+tkr 



(12) 



|k|<fc c 



Following Natoli and Ceperley [16], we express W as a linear combination of 
fifth-order polynomials, W(|r|) = Y^ n a nfn(\ r \), which is known as a locally 
piecewise-quintic Hermite interpolant. The fit coefficients a n and y k can be 
derived by minimizing the x 2 deviation, 



X 



dh 



V P (r)-J2a n f n (\r\)- ^ y k e + ^ 

n \k\<k c 



(13) 



where Q = L 3 is the volume of the unit cell. For the Fourier coefficients this 
directly yields, 



Vk = ^k - 5E a «/«k , with v k = ^ J dhV P (r)e lkr 



(14) 



where i>k and /„k are the corresponding Fourier transforms. Deviating from [16], 
we derive the coefficients a n from the following set of linear equations, m = 



6 



{l,...,n}, 



t'k/mk 

lkl<fe c 



fnm ^ ' fnkfmk 

|k|<fe c 



(15) 



v rn and / nm are overlap integrals, f nm = ± J Q d 3 r/„(r)/ m (r), and 

i / rf 3 r Vp(r) / m (|r|) = ± £ / rf 3 r W(r) / m (|r|) (16) 



n L n 



= ^Jdrr 2 V Y (r)f m (r) + J2^j drr fm(r) J dqqVy(q) . (17) 

L-r 

In the last expression, one must sum over a sufficiently large number of im- 
ages until the interaction is completely screened, |L| ^> D s . Computing the 
coefficients a n using the real-space integration in Eq. 17 is more efficient and 
accurate than the Fourier integration employed in [16]. Our approach also 
works well for the Coulomb problem, which was the motivation for the [16] 
work. In this case, the Ewald potential replaces V P in Eq. 13. 

This optimized Ewald approach provides us with an accurate representation 
of the periodic functions leading to efficient many-body simulations. We apply 
it to the Yukawa potential, Vy, to the corresponding pair action, Uy, and also 
to the kinetic energy term, (j^jf — Vy^j, unless it is very small for r = L/2. 
Typically, we use between 10 and 20 shells of k vectors. 



3 Results 



The phase diagram in Fig. 2 relates our simulations at a fixed density of r s = 
200 to the classical Yukawa melting computed by Hamaguchi, the classical 
Coulomb melting line given by T = - e k ^ T = 172, and the quantum melting for 
Coulomb systems [7]. 

The most straightforward way to detect melting in the simulation is to monitor 



the instantaneous value of the Lindemann ratio, 7 = y (u 2 )/rNN, which relates 
the average displacement of a particle from its original lattice site to the 
nearest neighbor distance. Fig. 3 shows the Lindemann ratio for three Monte 
Carlo simulations. At the beginning of each simulation, the particles are in the 
classical b.c.c. ground state. For temperatures sufficiently above the melting 
line, the system melts instantly. For temperatures only slightly above the 
melting line, the simulation shows a meta-stable superheated solid, which 
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Fig. 2. Phase diagram for Wigner crystal. The dashed lines show the classical melt- 
ing line, r = 172, and quantum melting line for Coulomb systems computed by 
Jones and Ceperley [7] with PIMC. The solid lines show the classical melting line 
for the Yukawa system with the screening length chosen for hydrogen. Our PIMC 
simulations at r s = 200 (o and o) show a significant reduction in the melting tem- 
perature below the classical value due to the quantum effects of the protons. 
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Fig. 3. The evolution of the Lindemann ratio is shown for the three Monte Carlo 
simulations in the vicinity of the melting temperature. 

might melt at some point during the MC simulation, as the black dashed line 
indicates. The time it takes to melt depends not only on temperature, but also 
on system size, the type of MC moves, and the MC random numbers, thereby 
making this criterion impractical for determining the melting temperature. 
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Fig. 4. Pair correlation functions, g(r), for different temperatures. 

Fig. 4 shows a series of pair correlation functions, g(r), for simulations at 
different temperatures. The magnitude of the oscillations in the g(r) show 
a significant temperature dependence in the liquid phase. However, there is 
only a small change upon melting and all g(r) functions in the solid phase are 
practically identical to the example shown for T=2.3xl0~ 5 . 
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Fig. 5. The sharp peaks in the structure factor, S(k), disappear when the system 
melts and a pattern typical for a liquid appears. This transition coincides with the 
increase in the Lindemann ratio. 

Compared to the pair correlation function, the structure factor, S(k), shows 
significant changes upon melting (Fig. 5). The disappearance of the peaks 
coincides with the increase in the Lindemann ratio beyond stable value of 
approximately 0.28. However, neither method can determine whether a sim- 
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ulation is in the meta-stable state. They only lead to an upper bound of the 
melting temperature. To determine the thermodynamic phase boundary, one 
needs the free energy in both phases, which can be obtained through thermo- 
dynamic integration of the internal energies. 
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T (Ha) 



Fig. 6. The internal energy per particle is shown as a function of temperature. The 
classical kinetic energy and the Madelung term have been removed. 

Fig. 6 compares PIMC internal energies with results from corresponding clas- 
sical MC that we have performed. At high temperature when the thermal de 
Broglie wavelength is short compared to the inter-particle spacing, the protons 
behave classically and the PIMC energies approach results from classical MC 
simulations. At low T, both result differ substantially due to the zero point 
motion. The analytical high T limit [14] is not reached exactly, due to finite 
size effects. 

In Fig. 7, we compare MC results with our results derived from the harmonic 
lattice approximation and the particle-in-a-cell (PIC) model. In the PIC ap- 
proximation, one assumes that the thermal motion of the particles are uncor- 
related. One freezes all particles in the supercell except one and derives all 
thermodynamic variables from the motion of this single classical particle. The 
PIC internal energies agree remarkably well with the corresponding classical 
MC results. A generalization to the quantum case is not straightforward. If 
one simply considers a quantum particle, represented by a path in a lattice 
of frozen particles then the resulting kinetic energies are far too high (worse 
than the harmonic approximation) because the remaining classical particles 
provide too high of a confining force due to the missing quantum fluctuations. 

Fig. 7 also shows results from the harmonic lattice approximation for a unit 
cell of corresponding size. Harmonic internal energies are significantly overesti- 
mated, primarily due to errors in the kinetic energy as demonstrated in Fig. 8. 
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T (Ha) 

Fig. 7. Low temperature region of Fig. 6. Results from the harmonic lattice approx- 
imation and from a classical particle-in-a-cell model have been added. Furthermore, 
PIMC and harmonic results for particles with mass=10 instead of 1 have been in- 
cluded to demonstrate that the harmonic approximation becomes more accurate as 
zero point fluctuations are reduced. 



0.1 




Fig. 8. Difference in internal, kinetic, and potential energy between PIMC and the 
harmonic lattice approximation. The sharp kinks near T=2.5xl0 -5 indicate the 
melting transition. 

The zero point motion of the protons is large enough so that paths travel into 
regions of the potential where the harmonic approximation is not longer valid. 
A comparison of the PIMC Lindemann ratios and the harmonic values shows 
that the harmonic approximation localizes the particles too much, thereby 
increasing the kinetic energy. 
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To further support this conclusion, we perform PIMC and harmonic calcula- 
tion with particles 10 times as heavy. Fig. 7 shows that the agreement improves 
substantially. 




Temperature T (Ha) 



Fig. 9. Internal and free energies, E and F, are shown for the solid and liquid phase. 

Fig. 9 shows the internal PIMC energies along with the free energies obtained 
from thermodynamic integration. The solid free energies agree very well with 
solid internal energies until melting, which suggest that each phonon mode 
is in its ground state, and excitation in the phonon spectra leads directly to 
melting. At low T, the free energies of the crystalline state are lower than the 
extrapolated values for the liquid, which means the solid phase is stable and 
the density of r s = 200 is not yet high enough to reach the quantum melting 
transition (see Fig. 2). 

At T=2.0xl0 -5 the free energies of both phases match, which determines the 
thermodynamic melting point. This corresponds to a temperature of 11,700 
K and a hydrogen mass density of 2,100 gem -3 . 

Fig. 9 also shows that a number of simulations that appeared to be stable were 
actually meta-stable, which can also occur in classical simulations (Fig. 7). 
The melting temperature of 2.0 xlO -5 is significantly below the corresponding 
classical value of 2.7x 10~ 5 , which suggests that quantum effects are important. 
However, a finite size extrapolation remains to be done. 



4 Conclusions 

In this article, we used many-body computer simulations of protons interacting 
via a Yukawa potential to model dense atomic hydrogen in the regime of 
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the Wigner crystal. Path integral Monte Carlo simulations were employed to 
capture the quantum effects of the protons. Electronic screening effects were 
treated in the Thomas-Fermi approximation, which distinguishes our results 
from the earlier work by Jones and Ceperley [7]. 

We use the Lindemann ratio, pair correlation functions, and the structure 
factor to study the stability of the Wigner crystal and to detect melting. We 
observed that the system can remain in a meta-stable state of a super-heated 
crystal during the entire course of a PIMC simulation, which makes a direct 
determination of the melting temperature very difficult. 

Instead, a reliable melting temperature can be obtained by matching of the 
free energies of both phases, which were derived by thermodynamic integration 
of the PIMC internal energies. For the density under consideration, r s = 200, 
we found that the quantum Yukawa systems melts at significantly lower tem- 
peratures than the corresponding classical system. 

Furthermore, we compared our PIMC results with other more approximate 
techniques. The harmonic lattice approximation overestimates the kinetic en- 
ergies significantly because the zero point motion of the protons is strong 
enough so that anharmonic effects in the crystal field become relevant. We also 
compared with a classical particle-in-a-cell model and found good agreement 
with classical MC simulation. However this method cannot be generalized 
easily to the case of quantum protons. 

We plan to extend our analysis to other densities and to derive a phase diagram 
that indicates the stability of the Wigner crystal of nuclei in the presence of 
electronic screening effects. A careful analysis of finite size effects also remains 
to be done. Future theoretical work on dense atomic hydrogen will need to 
describe the electronic properties on a more fundamental level. Coupled ion- 
electron Monte Carlo is one promising approach [17]. 
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